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We explore the complex dynamical behavior of simple predator-prey models of biological coevo- 
lution that account for interspecific and intraspecific competition for resources, as well as adaptive 
foraging behavior. In long kinetic Monte Carlo simulations of these models we find quite robust 
l//-like noise in species diversity and population sizes, as well as power-law distributions for the life- 
times of individual species and the durations of quiet periods of relative evolutionary stasis. In one 
model, based on the Holling Type II functional response, adaptive foraging produces a metastable 
low-diversity phase and a stable high-diversity phase. 
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I. INTRODUCTION 

Biological evolution presents a rich array of phenomena that involve nonlinear interactions between large numbers 
of units. As a consequence, problems in evolutionary biology have recently enjoyed increasing popularity among 
statistical and computational physicists. However, many of the models used by physicists have unrealistic features 
that prevent the results from attracting significant attention from biologists. In this paper we therefore develop and 
explore individual-based models of coevolution in predator-prey systems based on more realistic population dynamics 
than some earlier models. 0, H H H H HH H 

II. MODELS 

Recently the author, together with R. K. P. Zia, introduced a simplified form of the tangled-nature model 
of biological macroevolution, which was developed by Jensen and collaborators. H || In these simplified 
models, ja, @, 0, 0, llCt 1 1 1| the reproduction rates in an individual-based population dynamics with nonoverlap- 
ping generations provide the mechanism for selection between several interacting species. New species are introduced 
into the communit y th rough point mutations in a haploid, binary "genome" of length L, as in Eigen's model for molec- 
ular evolution. [E| [l3] The potential species are identified by the index / 6 [0, 2 L — 1]. (Typically, only N{t) <C 2 L 
of these species are present in the community at any one time t.) At the end of each generation, each individual of 
species / gives birth to a fixed number F of offspring with probability Pi before dying, or dies without offspring with 
probability (1 — Pi). Each offspring may mutate into a different species - generally with different properties - with a 
small probability fi. Mutation consists in flipping a randomly chosen bit in the genome. 

A. Simplified tangled-nature models 

In these models, the reproduction probability for an individual of species / is given by the nonlinear function 



where R is an external resource that is renewed at the same level each generation, and {nj(t)} is the set of population 
sizes of all the species resident in the community in generation t. The function A/ is given by 

A I (R, {nj(t)» = -h + Vl R/N tot {t) + M u nj(t)/N tot (t) - N tot (t)/N Q . (2) 

J 

Here bi is an "energy cost" of reproduction (always positive), and rji (positive for primary producers or autotrophs, and 
zero for consumers or heterotrophs) is the ability of individuals of sp ecies I to utilize the external resource R, while Nq 
is an environmental carrying capacity [l^j (a.k.a. Verhulst factor [Hf). The total population size is N to t(t) = ^ j nj(t). 
The main feature of this reproduction probability is the random interaction matrix M,|l6j which is constructed at 
the beginning of a simulation run, and thereafter kept constant (quenched randomness). If Mu is positive and Mji 
negative, / is a predator and J its prey, and vice versa. If both matrix elements are positive, the relationship is a 
mutualistic one, while both negative indicate an antagonistic relationship. 

Two versions of this model were studied in earlier work. In the first version, which we have called Model A, there is 
no external resource or birth cost, and the off-diagonal elements of M are stochastically independent and uniformly 
distributed over [—1,-1-1], while the diagonal elements are zero. This model evolves toward mutualistic communities, 
in which all species are connected by mutually positive interactions. 0, El EJ 

Of greater biological interest is a predator-prey version of the model, called Model B. In this case a small minority 
of the potential species (typically 5%) are primary producers, while the rest are consumers. The off-diagonal part 
of the interaction matrix is antisymmetric, with the additional restriction that a producer cannot also prey on a 
consumer. |8l 1 1 1| In simulations we have taken bi and the nonzero r\i as independent and uniformly distributed on 
(0, +1]. This model generates simple food webs with up to three trophic levels. pi UlL Xv\ 

Both of these models provide interesting results, which include intermittent dynamics with power spectral densities 
(PSDs) of diversities and population sizes that exhibit l//-like noise, as well as power-law distributions for the 
lifetimes of individual species and the duration of quiet periods of relative evolutionary stasis. From a theoretical 
point of view they also have the great advantage that the mean-field equation for the steady-state average population 
sizes in the absence of mutations reduces to a set of linear (if Nq = oo) or at most quadratic equations and thus can 
easily be solved exactly. [M ItI ITTl The models thus provide useful benchmarks for more realistic, but generally 
highly nonlinear models. 

In fact, the population dynamics defined by Eqs. Q and are not very realistic. In particular, by summing over 
positive and negative terms in A/, the models enable species with little food to remain near a steady state if they 
are also not very popular as prey, or have very low birth cost. A more serious problem is the ad-hoc nature of the 
normalization by the total population size N to t(t) in the resource and interaction terms in A/. While this is the source 
of the models' analytic solvability, it implies an indiscriminate, universal competition without regard to whether or 
not two species directly utilize the same resources or share a common predator. The purpose of the present paper is 
to develop models with more realistic population dynamics and explore the properties of their evolutionary dynamics. 



B. Functional-response model 



Here we develop a model with more realistic population dynamics that include competition between different 
predators that prey on the same species, as well as a saturation effect expected to occur for a predator with abundant 
prey. In doing so, we retain from the models discussed above the important role of the interaction matrix M, as well 
as the restriction to dynamics with nonoverlapping generations. 

We first deal with the competition between predator species by defining the number of individuals of J that are 
available as prey for /, corrected for competition from other predator species, as 

niM u 

where ^P rcd ( J ) runs over all L such that Mlj > 0, i.e., over all predators of J. Thus, ^P red ( J ) niJ = njj anc | jf / j s 
the only predator consuming J, then hu = nj. 

Analogously, we define the competition-adjusted external resources available to a producer species I as 

Rr = ^-R . (4) 
As in the case of predators, J^i Ri = R> an d a sole producer species has all of the external resources available to it: 
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Ri = R. With these definitions, the total, competition-adjusted resources available for the sustenance of species / are 



prcy(Z) 

Si = rjiRi + ^ Mijh u , (5) 
j 

where ^P rey ^ 7 - ) runs over all J such that Mjj > 0, i.e., over all prey of /. 

A central concept of the model is the functional response of species / with respect to J, <&/,/. 0,^1 This is the 
rate at which an individual of species / consumes individuals of J. The simplest functional response corresponds to 
the Lotka-Volterra model: [14( <&/j = nj if Mjj > and otherwise. However, it is reasonable to expect that the 
consumption rate should saturate in the presence of very abundant prey. |19| For ecosystems consisting of a single pair 
of predator and prey, or a simple chain reaching from a bottom-level producer through intermediate species to a top 
predator, the most common forms of functional response are due to Holling. 19] For more com plic ated, interconnected 
food webs, a number of functional forms have been proposed in the recent literat ureJ lSl 12(1 l2U I22L |22j but there is 
as yet no agreement about a standard form. Here we choose a ratio-dependent Holling Type II form.[T9l] 

Mijhu 

®u = . , (6) 

Abi + rij 

where A £ (0, 1] is the metabolic efficiency of converting prey biomass to predator offspring. Analogously, the functional 
response of a producer species I toward the external resource R is 

*IR = - • (7) 

ao j + nj 

In both cases, if ASi -C nj, then the consumption rate equals the resource (Mijhu or ijiRi) divided by the number 
of individuals of /, thus expressing intraspecific competition for scarce resources. In the opposite limit, XSi ^> nj, the 
consumption rate is proportional to the ratio of the specific, competition-adjusted resource to the competition-adjusted 
total available sustenance, Si. The total consumption rate for an individual of / is therefore 

prcy(J) ~ f ri t 

Cl -* IR+ \ ~ Js7Tn~i ~ I 1/A forA5 / »n 7 " (8) 
The birth probability is assumed to be proportional to the consumption rate, 

Bj = A<7je[0,+l], (9) 
while the probability that an individual of / avoids death by predation until attempting to reproduce is 

prcd(J) 

Ai = l- V (10) 
The total reproduction probability for an individual of species I in this model is thus Pi it) — Aiit)Biit). 



III. NUMERICAL RESULTS FOR THE FUNCTIONAL-RESPONSE MODEL 

We simulated the functional-response model over 2 24 = 16 7 7 7 2 1 6 generations (plus 2 20 generations "warm-up") 
for the following parameters: genome length L = 21 (2 21 = 2 097 152 potential species), external resource R — 16 000, 
fecundity F = 2, mutation rate [i — 10~ 3 , proportion of producers c pro d = 0.05, interaction matrix M with connectance 
C = 0.1 and nonzero elements with a symmetric, triangular distribution over [— 1,+1], and A = 1.0. We ran five 
independent runs, each starting from 100 individuals of a single, randomly chosen producer species. 



A. Time series 



Time series of diversities (effective numbers of species) and population sizes for one realization are shown in Fig.^ To 
filter out noise from low-population, unsuccessful mutations, we define the diversity as the exponential Shannon- Wiener 
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FIG. 1: (Color online.) (a) Time series of diversities (measured in species, lower curves) and population sizes (measured in 
individuals, upper curves) for one specific simulation run. The strongly fluctuating curves in the background are sampled every 
8192 generations, while the smooth curves with data points in contrasting colors that are overlaid in the foreground are running 
averages over 524 288 generations. Black with light gray (yellow online) overlay: all species. Light gray with dark gray overlay 
(green and magenta online): producers. Dark gray with light gray overlay (red and cyan online): consumers, (b) Time series of 
the accumulation of new and extinct species in the same simulation run depicted in (a) . The solid, black curve shows the total 
number of different species that have at least once attained a population size ni > 1000 by time t. The dashed curves count 
the total number of species that have gone extinct after attaining a maximum population greater than 1000. The black dashed 
curve refers to all species, the light gray one (green online) to producers, and the dark gray one (red online) to consumers. The 
ratio of approximately 1.89 between the dashed and full black curves indicate that major species recur on average about twice 
during the evolution. This is an artifact of the finite genome length. The inset shows the detailed, intermittent structure of 
the solid, black curve over 400 000 generations. The interval is indicated by a horizontal bar in the main panel. 



index. 26] This is the exponential function of the information-theoretical entropy of the population distributions, 
D(t) = exp [S ({ni(t)})], where 

S({m(t)}) = - ]T Pi(t)ln Pl (t) (11) 

{I\pi(t)>0} 

with pi(t) = ni(t)/N to t(t) for the case of all species, and analogously for the producers and consumers separately. 

The time series for both diversities and population sizes show intermittent behavior with quiet periods of varying 
lengths, separated by periods of high evolutionary activity. In this respect, the results are similar to those seen for 
Models A and B in earlier work. |ol M IsL fTol ITl| However, diverse communities in this model seem to be less stable 
than those produced by the linear models. In particular, this model has a tendency to flip randomly between an active 
phase with a diversity near ten, and a "garden of Eden" phase of one or a few producers with a very low population 
of unstable consumers, such as the one seen around 10 million generations in Fig. ^ 

B. Power-spectral densities 

A common method to obtain information about the intensity of fluctuations in a time series at different time scales 
is the power-spectral density (squared Fourier transform), or PSD. PSDs are presented in Fig. [21 for the diversity 
fluctuations and the fluctuations in the population sizes (Fig.[3{a)) and the intensity of extinction events (Fig. Gib)). 
The former two are shown for the total population, as well as separately for the producers and consumers. All three 
are similar. Extinction events are recorded as the number of species that have attained a population size greater 
than one, which go extinct in generation t (marked as "species" in the figure), while extinction sizes are calculated 
by adding the maximum populations attained by all species that go extinct in generation t (marked as "population" 
in the figure). The PSDs for all the quantities shown exhibit approximate 1/f behavior. For the diversities and 
population sizes, this power law extends over more than five decades in time. The extinction measures, on the other 
hand, have a large background of white noise for frequencies above 10 -3 generations -1 , probably due to the high rate 
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FIG. 2: (Color online.) (a) PSDs for the diversities and population sizes, each recorded separately for all species and for 
producers and consumers. The time series were sampled every 16 generations, (b) PSDs for the extinction activity. In both 
parts of the figure, the results are averaged over five independent simulation runs. See discussion in the text. 
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FIG. 3: (Color online.) (a) Histograms of species lifetimes, shown for all species, as well as separately for producers and 
consumers, (b) Histograms of the durations of evolutionarily quiet periods, defined as the times that |dS/dt| (averaged over 
16 generations) falls continuously below some cutoff. The inset is a histogram of d5*/dt, showing a Gaussian center with 
approximately exponential wings. The parabola in the foreground is a Gaussian fit to this central peak. The cutoff values for 
the main figure, between 0.008 and 0.024, were chosen on the basis of this distribution. The data in both parts of the figure 
are averaged over five independent simulation runs. 



of extinction of unsuccessful mutants. For lower frequencies, however, the behavior is consistent with 1// noise within 
the limited accuracy of our results. 



C. Species lifetimes and durations of quiet periods 



The evolutionary dynamics can also be characterized by histograms of characteristic time intervals, such as the 
time from creation till extinction of a species (species lifetimes) or the time intervals during which some indicator of 
evolutionary activity remains continuously below a chosen cutoff (duration of evolutionarily quiet periods) . Histograms 
of species lifetimes are shown in Fig. 01a) . As our indicator of evolutionary activity we use the magnitude of the 
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FIG. 4: (Color online.) (a) Time series of population sizes (upper curves) and diversities (lower curves) for the model with 
adaptive foraging. The interpretation of the colors and lines are the same as in Fig. Eta), (b) Time series of the number of new 
species that have reached a population greater than 1000 (lower curve) and greater than 100 (upper curve). The inset shows 
the intermittent structure of the upper curve on a very fine scale of 2000 generations. See discussion in the text. 



logarithmic derivative of the diversity, |dS/di|, and histograms for the resulting durations of quiet periods, calculated 
with different cutoffs, are shown in Fig. E^b). Both quantities display approximate power-law behavior with an 
exponent near —2, consistent with the 1// behavior observed in the PSDs.[f|, [21} It is interesting to note that the 
distributions for these two quantities for this model have approximately the same exponent. This is consistent with 
the previously studied, mutualistic Model A, 0,0 but not with the predator-prey Model B.p. lllllr^ | We believe the 
linking of the power laws for the species lifetimes and the duration of quiet periods indicate that the communities 
formed by the model are relatively fragile, so that all member species tend to go extinct together in a "mass extinction." 
In contrast, Model B produces simple food webs that are much more resilient against the loss of a few species, and as 
a result the distribution of quiet-period durations decays with an exponent near — l.pj ITU Il7| 



IV. ADAPTIVE FORAGING 



The model studied above is one in which species forage indiscriminately over all available resources, with the 
output only limited by competition. Also, there is an implication that an individual's total foraging effort increases 
proportionally with the number of species to which it is connected by a positive M/j. A more realistic picture 
would be that an individual's total foraging effort is constant and can either be divided equally, or concentrated on 
richer resources. This is known as adaptive foraging. While one can go to great length devising optimal foraging 
strategies. |l8l 12^ we here only use a simple scheme, in which individuals of / show a preference for prey species J, 
based on the interactions and population sizes (uncorrected for interspecific competition) and given by 

Mijnj 
thR + YJk m ikti k 



and analogously for R by 



niR 



9IR = „ , ^prey(/) " ( 13 ) 

wR+2Jk M IK n K 

The total foraging effort is thus gm + g TJ = \, The preference factors are used to modify the reproduction 

probabilities by replacing all occurrences of Mjj by Mjjgjj and of m by mgiR in Eqs. QHZ). 

The results of implementing the adaptive foraging are quite striking. The system appears now to have a metastable 
low-diversity phase similar to the active phase of the non-adaptive model, from which it switches at a random time 
to an apparently stable high-diversity phase with much smaller fluctuations. As seen in Fig.^Ja), the switchover is 
quite abrupt, and Fig.@£b) shows that it is accompanied by a sudden reduction in the rate of creation of new species. 
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FIG. 5: (Color online.) PSDs for the diversities and population sizes in the metastable phase (averaged over three independent 
runs) (a), and the stable phase (averaged over five independent runs) (b) for the model with adaptive foraging. Both show 
approximate 1/f noise for frequencies above about 10 , but the PSDs appear to approach constant levels for the lowest 
frequencies. 

As seen in Fig. [SJ the PSDs for both the diversities and population sizes in both phases show approximate 1/f noise 
for frequencies above 10 -5 generations -1 . For lower frequencies, the metastable phase shows no discernible frequency 
dependence, while for the stable phase, the frequency dependence continues at least another decade. It thus appears 
that long-time correlations are not seen beyond 10 generations for the metastable phase, and probably not beyond 
about 10 6 generations for the stable one. These observations are consistent with species-lifetime distributions for both 
phases (not shown), which are quite similar to those for the non- adaptive model, but typically with cutoffs in the 
range of 10 5 to 10 6 generations, much shorter than the total simulation times. 

In fact, the system can also escape from the low-diversity phase to total extinction, which is an absorbing state, 
and in some of our simulation runs we avoided this by limiting |Afjj| to less than 0.9. This restriction does not seem 
to have any effect on the dynamics in the high-diversity phase. These results are preliminary, and it is possible that 
the high-diversity phase corresponds to a mutational meltdown. More research is clearly needed regarding the effects 
of adaptive foraging in this model. 

V. CONCLUSIONS 

In this paper we have shown that very complex and diverse dynamical behavior results, even from highly over- 
simplified models of biological macroevolution. In particular, PSDs that show l//-like noise and power-law lifetime 
distributions for species as well as evolutionarily quiet states are generally seen. This is the case, both in the 
analytically tractable, but somewhat unrealistic tangled-nature type models, and in the nonlinear predator-prey 
models based on the more realistic Holling Type II functional response. Particularly intriguing is the appearance of a 
new, stable high-diversity phase in the latter type of model when adaptive foraging behavior is included. Among the 
many questions about this new phase that remain to be addressed is the structure of the resulting community food 
webs. 
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